SECONDARY  CIRCULATION  IN  GRANULAR  FLOW  THROUGH 
NONAXISYMMETRIC  HOPPERS 
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Abstract.  Jenike’s  radial  solution,  widely  used  in  the  design  of  materials-handling  equipment, 
is  a  similarity  solution  of  steady-state  continuum  equations  for  the  flow  under  gravity  of  granular 
material  through  an  infinite,  right-circular  cone.  In  this  paper  we  study  how  the  geometry  of  the 
hopper  influences  this  solution.  Using  perturbation  theory,  we  compute  a  first-order  correction  to 
the  (steady-state)  velocity  resulting  from  a  small  change  in  hopper  geometry,  either  distortion  of  the 
cross  section  or  tilting  away  from  vertical.  Unlike  for  the  Jenike  solution,  all  three  components  of  the 
correction  velocity  are  nonzero:  i.e.,  there  is  secondary  circulation  in  the  perturbed  flow.  We  show 
that,  depending  on  hopper  and  material  parameters,  the  perturbed  velocity  depends  sensitively,  to  an 
astonishing  degree,  on  hopper  geometry.  These  results  suggest  that,  even  in  a  vertical  conical  hopper, 
solutions  with  circulation  may  bifurcate  from  the  Jenike  solution,  a  phenomenon  to  be  investigated 
in  a  future  paper. 
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1.  Introduction.  In  manufacturing  industries,  raw  materials  are  stored  in  gran¬ 
ular  form  in  a  silo,  and  when  needed,  they  are  expected  to  flow  out  of  the  silo  under 
gravity  through  a  hopper.  Problems  in  the  discharge  process  are  frequent  and  expen¬ 
sive,  see  e.g.  [9].  As  demonstrated  by  a  Rand  Corporation  study  [10],  these  problems 
are  symptomatic  of  our  poor  understanding  of  the  behavior  of  granular  materials1. 

Jenike’s  radial  solution  is  a  central  component  of  silo  design.  Despite  its  impor¬ 
tance,  this  solution  is  subject  to  many  severe  restrictions: 

1.  Granular  material  is  modeled  as  a  continuum,  with  an  ad  hoc  constitutive  law. 

2.  The  flow  is  assumed  to  be  steady. 

3.  The  flow  domain,  a  mathematical  idealization,  is  an  infinite  cone,  given  in 
spherical  polar  coordinates  by  the  formula 

{(r,  0, 4>)  :  0  <  r  <  oo,  0  <  6  <  6W },  (i 9W  =  constant). 

4.  Only  similarity  solutions  are  considered. 
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1The  study  compared  the  design  output  and  the  actual  output  of  a  total  of  60  manufacturing 
plants  in  various  industries,  22  that  were  based  primarily  on  liquids-processing  technology  and  38  on 
solids-processing  technology.  On  average,  the  liquids-processing  plants  produced  at  84%  of  design 
capacity  while  the  solids-processing  plants  produced  at  only  63%  of  design  capacity.  To  quote  Merrow, 
“In  economic  terms,  the  difference  between  63%  of  design  and  84%  is  very  large.  It  implies  a  capital 
cost  per  unit  of  output  about  one-third  higher  for  the  solids-processing  plants,  on  the  basis  of  poor 
performance  alone.  In  addition,  poor  performance  is  inevitably  associated  with  higher  operating  and 
maintenance  costs  per  unit  of  product.”  Moreover,  the  standard  deviation  of  the  solids-processing 
data  set  was  much  greater,  indicative  of  our  difficulties  in  predicting  the  behavior  of  granular  solids. 
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In  this  paper  we  relax  restrictions  3  and  4  partly.  Specifically,  we  generalize  the 
domain  to  an  infinite  pyramidal  hopper  described  by  the  inequality 

(1.1)  0  <  6  <  9W  +  e  cos  mcj), 

where  e  is  a  small  parameter  and  m  is  a  positive  integer.  Assuming  a  perturbation 
series 


v<°>  +  +  . . . 

for  the  flow  velocity  in  the  domain  (1.1),  where  v is  Jenike’s  solution,  we  derive 
a  linear  PDE  for  the  first-order  correction  t/P.  The  r-dependence  of  i/P  still  has 
similarity  form,  and  the  ^-dependence  may  be  handled  by  separation  of  variables.  In 
this  way  we  reduce  solving  the  PDE  for  r/P  to  solving  a  two-point  boundary  problem 
on  the  interval  0  <  9  <  9W. 

In  Jenike’s  solution,  only  the  radial  component  Ur0'*  of  the  velocity  is  nonzero.  By 
contrast,  all  three  components  of  the  correction  velocity  t/P  are  nonzero.  In  other 
words,  distortion  of  the  conical  domain  leads  to  secondary  circulation.  For  example,  in 
Figure  5.1  below,  the  flow  in  the  6,  (/-directions  is  shown  in  two  cases  which  correspond 
to  a  circular  hopper  that  is  tilted  slightly  to  the  right,  and  in  Figure  5.2,  in  two  cases 
which  correspond  to  a  slightly  distorted  vertical  hopper. 

The  boundary  problem  for  i/P  contains  three  significant  parameters:  the  angle 
of  internal  friction  S,  the  coefficient  of  wall  friction  /i„, ,  and  the  opening  angle  of  the 
hopper  6W.  (The  subscript  w  is  mnemonic  for  wall.)  Surprising  behavior  occurs  when 
these  parameters  are  varied.  In  the  first  place,  the  direction  of  circulation  may  reverse 
itself.  This  is  illustrated  for  instance  by  Figure  5.2:  6  and  9W  are  the  same  for  both 
parts  of  the  figure,  but  pw  is  larger  in  the  bottom  part.  As  illustrated  in  Figure  5.1, 
the  topology  of  the  circulation  may  change  as  fiw  varies.  Most  surprising  of  all,  the 
circulation  does  not  reverse  direction  by  passing  smoothly  through  zero;  rather,  as 
illustrated  by  the  graph  of  v ^  in  Figure  5.3,  the  circulation  suffers  a  “1/x-blowup” 
as  the  parameters  pass  through  the  critical  values!  (Of  course  in  deriving  the  PDE  for 
i/P,  it  was  assumed  that  r/P  was  much  smaller  than  v^°\  When  this  PDE  predicts 
that  i/P  is  large,  the  derivation  fails.  Thus,  the  blowup  of  r/P  does  not  mean  that  the 
solution  of  the  full  nonlinear  problem  diverges,  but  it  does  mean  that  the  full  solution 
does  not  depend  smoothly  on  parameters.  In  other  words,  the  nonlinear  solution  may 
be  extremely  sensitive  to  perturbations.) 

The  outline  of  the  paper  is  as  follows.  In  Section  2,  the  governing  equations 
are  recalled  together  with  Jenike’s  construction  of  similarity  solutions  in  conical  do¬ 
mains.  For  nonaxisymmetric  domains  of  the  type  (1.1),  the  problem  is  then  linearized 
about  Jenike’s  solution  in  Section  3.  The  resulting  system  is  discretized  in  Section  4. 
Numerical  results  and  discussion  are  offered  in  Section  5. 

2.  The  model. 

2.1.  Governing  equations  and  boundary  conditions.  The  unknowns  are 
the  3-component  velocity  vector  v,  the  3x3  symmetric  stress  tensor  T,  and  a  scalar 
plasticity  coefficient  A.  (The  density  p  is  a  constant.)  In  total,  there  are  3+6+1=10 
unknown  functions.  In  writing  the  equations  for  these  variables,  we  need  the  strain 
rate  tensor  V  =  —1/2 (Vu  + VuT)  and  the  deviatoric  part  of  the  stress  tensor  devT  = 
T  —  ItrTI.  Note  the  sign  convention:  V  measures  the  compression  rate  of  the 
material;  analogously,  positive  eigenvalues  of  T  correspond  to  compressive  stresses. 
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This  sign  convention  reflects  the  fact  that  granular  materials  disintegrate  under  tensile 
stresses. 

Following  [12],  we  require  that  these  variables  satisfy 

(2.1)  V-T  =  p5, 

(2.2)  V  =  X  dev  T, 

(2.3)  |devT|2  =  2s2(trT/3)2, 

where  g  is  the  (vector)  acceleration  of  gravity,  |  •  |  denotes  the  Frobenius  norm 

iti2=  E  Tfj  =  ^T2 

i,j=l 

(the  latter  equality  only  for  symmetric  tensors)  and  s  =  sin  6,  with  <5  being  the  angle 
of  internal  friction  of  the  material  under  consideration  (see  [11]).  Equation  (2.1) 
expresses  force  balance:  i.e.,  Newton’s  second  law  with  inertia  neglected  because  the 
flow  is  assumed  slow;  it  is  equivalent  to  three  scalar  equations.  Equations  (2.2)  and 

(2.3)  are  constitutive  laws,  the  alignment  condition  and  the  von  Mises  yield  condition, 
respectively;  they  are  equivalent  to  six  and  to  one  scalar  equations,  respectively.  Thus 
(2. 1-2.3)  is  a  determined  system,  10  equations  for  10  unknowns.  Since  (2.3)  contains 
no  derivatives,  this  system  has  a  differential-algebraic  character.  Taking  the  trace  of 
(2.2),  we  see  that  div  v  =  — tr  V  =  0;  thus,  incompressibility  is  part  of  the  constitutive 
assumptions.  Incidentally,  for  a  solution  to  be  physical,  the  function  A  in  (2.2)  must 
satisfy  A  >  0  everywhere;  otherwise  friction  would  be  adding  energy  to  the  system 
rather  than  dissipating  it.  In  fact,  we  want  A  to  be  strictly  positive  since  one  of 
the  assumptions  underlying  the  derivation  of  (2. 1-2. 3)  is  that  material  is  actually 
deforming. 

We  seek  solutions  of  (2. 1-2.3)  in  a  pyramidal  domain,  expressed  in  spherical  polar 
coordinates  as 

(2.4)  n  =  {(r,0,0):O  <6<C(<f>)}, 

where  C  is  a  given  smooth  27r-periodic  function.  Such  a  domain  represents  a  mathe¬ 
matical  idealization  of  a  converging  hopper,  in  general  a  nonaxisymmetric  one. 

On  the  boundary  dfl  =  {(r,C(<j)),  </>)},  wall  impenetrability  imposes  one  boundary 
condition  on  the  velocity:  i.e., 

(2.5)  vN  =  0, 

where  Vn  is  the  normal  velocity.  Two  additional  boundary  conditions  come  from 
Coulomb’s  law  of  sliding  friction.  The  surface  traction  r — i.e.,  the  force  exerted  by 
the  wall  on  the  material — is  given  by 

3 

Ti  =  E  Nj , 

3  =  1 

where  N  is  the  unit  interior  normal  to  d Q.  If  the  vector  r  has  normal  component  tn 
and  tangential  component  tt  =  t  —  t^N,  then  we  require  that 


(2.6) 


tt  =  -HwTN{v/\v\) 
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where  /i,„  is  the  coefficient  of  friction  between  the  wall  and  the  material.  Note  that: 
(i)  If  T  is  positive  definite  (i.e.,  if  all  stresses  are  compressive),  then  tn  >  0.  (ii)  While 
tn  is  a  scalar,  tt  is  effectively  a  two-component  vector;  thus,  (2.6)  is  equivalent  to 
two  scalar  equations,  (iii)  Because  of  (2.5),  the  velocity  v  is  tangential  to  dih  we  are 
assuming  that  v  ^  0  at  the  boundary. 

2.2.  Jenike’s  similarity  solution.  Suppose  that  the  domain  (2.4)  is  axisym- 
metric:  i.e.,  suppose 

(2.7)  n  =  {(r,6,(f>)  -0<0<9W}, 

where  9W  is  a  constant.  In  this  case  Jenike  [8]  found  that  (2. 1-2. 3)  have  solutions  that 
are  independent  of  </>  and  have  a  similarity  dependence  on  r, 

v(o)  e)=r~2  t)(0)  (0) ,  T(0)  (r,  9)  =  r  f (0)  (0) . 

(Here  and  below,  a  hat  above  a  variable  indicates  a  function  that  depends  on  9  alone.) 
Moreover,  only  the  radial  component  of  velocity  is  nonzero:  i.e.,  Vg =  0. 
Similarly  =  Tg^  =  0.  Indeed  all  components  of  T  can  be  expressed  in  terms  of 
two  scalar  variables,  the  so-called  Sokolovskii  variables  [11],  the  mean  stress  p ^  = 
tr  T ^  /3  and  an  angle  ip'i  specifically, 

—  cos  2 ip  —  sin  2ip  0 

(2.8)  devT(0)  =  sp(0)  -sin2^>  ^cos2  0  , 

0  0  cos  2 ip 

where  p ®  =  rp^  and  the  function  if),  like  p^,  depends  only  on  9. 

The  boundary  conditions  (2. 5, 2. 6)  may  be  written  more  explicitly  when  f 1  is 
axisymmetric.  Equation  (2.5)  reduces  to 

(2.9)  ve  =  0. 

Let  us  decompose  the  vector  equation  (2.6)  into  a  direction  and  a  magnitude.  Re¬ 
garding  the  direction,  the  vectors  tt  and  v  are  parallel  if 

(2.10)  Tr^Vf/,  -  Tg^Vr  =  0. 

Jenike’s  solution  satisfies  both  (2.9)  and  (2.10)  trivially.  The  two  sides  of  (2.6)  have 
equal  magnitude  if 

(2.11)  Trg  +  pwTgg  =  0. 

We  briefly  summarize  the  construction  of  Jenike’s  solution,  referring  to  [12]  for 
more  details.  The  ansatz  (2.8)  arranges  that  (2.3)  holds  automatically.  On  substitu¬ 
tion  into  (2.1),  we  obtain  a  first-order  2x2  system  of  ordinary  differential  equations 
for  p*'0)  and  if '>.  This  system  has  a  regular  singular  point  at  0  =  0,  and  one  boundary 
condition  comes  from  requiring  that  the  solution  be  regular  there;  the  other  boundary 
condition  comes  from  (2.11).  Thus,  the  stresses  are  determined  as  the  solution  of  a 
two-point  boundary-value  problem.  (In  axial  symmetry,  the  stress  equations  decouple 
from  the  velocity.)  Once  the  stresses  are  known,  (2.2)  reduces  to  a  linear  first-order 
ODE  for  i>r.  The  velocity  is  determined  only  up  to  a  multiplicative  constant,  but 
the  normalization  of  the  velocity  will  scale  out  of  the  calculations  below. 
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Incidentally,  for  Jenike’s  solution  the  plasticity  coefficient  A  in  (2.2),  which  cancels 
out  in  the  derivation  of  the  equation  for  vl0^ ,  has  the  form 

A(0)(r,6»)  =r"4A(0)(6»). 

Using  (2.2),  the  function  A*'0-*  may  be  determined  from  vl°\ 

3.  Linearized  analysis  for  a  nearly  axisymmetric  domain. 

3.1.  Derivation  of  linearized  differential  equations.  Equations  (2. 1-2. 3), 
a  10  x  10  nonlinear  DAE  system  that  is  elliptic  in  the  sense  of  Agmon,  Douglis,  and 
Nirenberg  [1],  present  formidable  mathematical  and  numerical  challenges.  In  this 
paper,  we  consider  a  simplified  problem  that  exhibits  some  astonishing  behavior  of, 
and  prepares  the  way  for  computations  with,  the  full  problem  on  a  general  domain. 
Suppose  the  function  C  specifying  the  boundary  of  fl  in  (2.4)  has  the  expansion 

(3.1)  C  (</))=  9W  +  e  cos  (to  <j>)  +  0(e2) 

where  m  is  a  positive  integer.  For  example,  a  slightly  tilted  (circular)  cone  admits 
such  a  representation  with  m  =  1,  where  e  measures  the  angle  of  tilt;  likewise  for  a 
(vertical)  pyramidal  hopper  having  a  slightly  elliptical  cross  section,  with  m  =  2. 

An  expansion  of  the  solution 

(3.2)  u  =  u(0)+eu(1)  +  C>(e2),  T  =  T(0)  +  eT«  +  0(e2) 

is  sought,  where  i are  equal  to  Jenike’s  radial  solution  [8].  Substituting  (3.2) 
into  (2. 1-2. 3),  we  derive  the  equations  for  the  first-order  perturbation 

(3.3)  V  •  T(1)  =  0, 

(3.4)  V(1)  =  A(1)  devT(0)  +  A(0)devT(1) , 

(3.5)  tr(devT(0)  devT(1))  =  2  s2  p(0)p(1) 

where  =  tr  TW /3,  i  =  0, 1  are  the  mean  stresses. 

The  correction  velocity  v W  has  the  same  r-dependence  as  the  Jenike  solution 
[8]  (although  all  three  components  of  v ^  are  nonzero),  and  its  ^-dependence  can  be 
obtained  through  separation  of  variables.  Indeed,  suppose  each  component  of  has 
the  form 

(3.6)  v!j ^  =  r~2  ij ^  ( 9 )  trig  mcj) 

where  trig  m<p  denotes  either  cos  mxj>  or  sin  mxj>.  In  order  to  satisfy  the  appropriately 
modified  version  of  the  boundary  condition  (2.9)  on  the  perturbed  domain,  v ^  will 
have  to  be  in  phase  with  (3.1):  i.e.,  we  need 

=  r~2  6^  ( 9 )  cos  m.cj). 

It  is  readily  seen  that  if 

=  r~2  ( 9 )  cos  mtj)  and  =  r~2  ( 9 )  sin  m<t>, 

then  all  terms  in 

(3.7)  V  •  =  drv^  +  2r~1v^  +  r~1d$v ^  +  r^1  cot  Ov fj1^  +  r_1  esc  dd^v^ 
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Scalars:  p  =  rp(9)  cos  rruf) 


Vectors:  v  =  \ 


vr{6)  cos  rruf) 
vg(0)  cos  rruf) 
v<l>(9)  sin rruf) 


Tensors: 


T  =  r 


Trr(9)  cos  rruf) 
Trg  ( 9 )  cos  rruf) 
Tgg  (9)  cos  rruf) 
Tr${9)  sin  rruf) 
Tg</,(9)  sin  rruf) 
T^(9)  cos  rruf) 


Table  3.1 

The  r-  and  ^-dependence  of  scalars,  vectors,  and  tensors  in  separation  of  variables 


are  proportional  to  r  3  cos  i.e.,  variables  separate  in  the  equation  V  •  v  =  0. 

Tables  3. 1-3.3  help  systematize  the  elimination  of  ^-dependence  in  (3. 3-3. 5)  with 
separation  of  variables.  The  appropriate  r-  and  (^-dependence  for  the  scalar  p^\  for 
the  vector  v^\  and  for  the  tensor  T ^  is  indicated  in  Table  3.1.  (Note  that  symmetric 
3x3  tensors  are  represented  as  vectors  in  R6,  the  components  being  enumerated  in 
the  order  shown.)  In  Table  3.2  we  record,  for  the  reader’s  convenience,  the  expressions 
in  spherical  coordinates  for  four  differential  operators  that  occur  in  these  equations. 


Vp=[drp,  r  1dgp,  r  1  esc.  9d<t,p]T 


V  ■  v  =  drvr  +  2r  1vr  +  r  1dgvg  +  r  1  cot  9vg  +  r  1  esc 


br 

dr  Vf 

VrO 

\  ( r~1dgvr  -  r~xvg  +  drVg) 

Vee 

r-1  ( Vr  +  dgvg) 

Vrrfi 

\  (r-1  esc  9  d,j,Vr  —  r~1vcf,  +  drv $) 

Ve<t> 

|r_1  (dgi v,  -  cot  9  v#  +  esc  9d^vg) 

r_1  (vr  +  cot  9vg  +  esc  9  d^v^) 

V  •  T  = 


drTrr  +  r  1  esc  9  d<f,Tr4>  +  r  1  dgTrg  +  r  1  (2Trr  —  —  Tgg  +  Trg  cot  9 ) 

drTre  +  r_1  esc  9  d^Tg^  +  r~1dgTgg  +  r”1  (3 Trg  +  ( Tgg  —  cot  9) 
drTr(j,  +  r-1  csc9  +  r~1dgTg<j>  +  r_1(3Tr^  +  2Tg^,  cot  6) 


Table  3.2 

Differential  operators  in  spherical  coordinates 


The  main  point,  which  makes  separation  of  variables  work  in  this  problem,  is  that 
the  0-dependent  part  of  each  of  these  linear  operators  is  given  by 


(3.8) 


(a) 

Vp 

=  ( gide  +  go)p 

(b) 

V  •  v 

=  (djde  +  d^)v 

(c) 

V 

=  -(Gtfe  +  Go)* 

(d) 

VT 

=  (D1dg  +  D0)T 

where  gi,  go, . . . ,  D0  are  the  matrices  given  in  Table  3.3. 

The  calculation  needed  to  verify  (3.8b)  was  described  above;  the  other  equations 
may  be  verified  similarly.  Incidentally,  (3.8a)  may  be  derived  by  substituting  T  =  pi 
in  (3.8d),  and  (3.8b)  may  be  derived  by  taking  the  trace  of  (3.8c). 
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9i  = 

=  [  0 

1 

0 

]T 

So  : 

=  [  1 

0 

sin  o 

r 

d  i  = 

=  [  0 

1 

0 

]T 

do  = 

[  0 

cote  1 T 

sin  9  J 

'  0 

0 

0  " 

r 

-2 

0 

0  " 

1/2 

0 

0 

0 

-3/2 

0 

G\ 

0 

1 

0 

Go 

1 

0 

0 

— 

0 

0 

0 

— 

- 

m 

2  sin  6 

0 

-3/2 

0 

0 

1/2 

0 

m 

2  sin  0 

cot  6 

2 

0 

0 

0 

l 

cot  9 

m 

sin  6 

'  0 

1 

0 

0 

0 

0  " 

'  3 

cot  6 

-l 

m 

sin  0 

0 

-1 

D1  = 

0 

0 

1 

0 

0 

0 

Do  = 

0 

4 

cot  6 

0 

m 

sin  0 

-  cot  8 

0 

0 

0 

0 

1 

0 

0 

0 

0 

4 

2  cot  6 

m 

sin  9 

Table  3.3 
Matrices  in  (3.8) 


With  this  notation,  (3. 3-3. 5)  reduces  to  a  system  of  ODEs  in  0, 

(3.9)  (Dtfe  +  D0)T «  =  0, 

(3.10)  -{Gide  +  G'o)0(1)  =  A(1)  devT^  +  A^devT^, 

(3.11)  tr(devT(0)devT(1))  =  2  s2p(0)p(1). 

Recalling  the  representation  of  symmetric  tensors  as  6-component  vectors,  we  observe 
that  the  LHS  of  (3.11)  may  be  rewritten  as  an  inner  product 

tr  (devT(0)devT(1))  =  devf(0)TAl  devf(1) 

where  A4  is  the  6x6  matrix 


M  =  diag(l,  2, 1,  2,  2, 1); 
thus,  we  may  rewrite  (3.11)  as 

(3.12)  clevf  (0)TAd  devf (1)  =2  s2p(0)p(1). 


Let  us  show  that  the  deviatoric  stresses  in  (3.9,3.10,3.12)  can  be  eliminated  from 
these  equations  to  obtain 


(3.13)  (B1de  +  D0)  ^-J-(Glds  +  G0)v(1)  -  ^^(0)  j  +  (. gidg  +  go)p{1)  =  0 

(3.14)  (dfae  +  4)t)(1)  =  0 


where  in  (3.13) 

A«  = 


2s2  (p(0))2^(0) 


V^T  M(Gidg  +  Go)i 


,(i)  _ 


AW 

»(°) 


di) 


(3.15) 
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Equation  (3.14)  follows  on  taking  the  trace  of  (3.10).  Next,  we  rewrite  (3.10)  as 

1  \(i) 

(3.16)  devfW  =  (Gido  +  G0)u(1)  -  C(0) 

v  '  AW  (A(°))2 

where  we  have  eliminated  devT^0)  using  the  relation  (7°)  =  A^devT^0) — effectively, 
equation  (2.2)  for  Jenike’s  solution.  Recalling  that  =  devT ^  +  p^I,  we  substi¬ 
tute  (3.16)  into  (3.9)  to  derive  (3.13).  Similarly,  (3.15)  follows  on  substituting  (3.16) 
into  (3.12)  and  rearranging. 

As  a  final  simplification,  we  substitute  (3.15)  into  (3.13),  obtaining  the  linear, 
homogeneous  system  of  ODEs 

(3.17)  -(A2dee  +  A\dg  +  A0)n(1)  +  (b^dg  +  b0)p(1)  =  0 

(3.18)  (df0j  +  =  0 


where,  with  the  definition 

P=  I  - 


1 


2s2(p(°))2(A(°))2 
the  coefficient  matrices  are  given  by 

A2  = 

A(°) 

Ai  =  jj^(D0PGi  +  DiPGq)  +  Didg 


v(°)\ X°)tm, 


VA(°) 


PG, 


^4°  —  W)D«PGo  +  Didg 


(J—PG( 

Va(°) 


bi  =  9i  +  Di 


T>(°) 

p(°)  A(°) 


bo  —  9o  +  {D\dg  +  D0) 


V(o) 


These  matrices  depend  on  9  and  in  fact  are  singular  as  9  — »  0.  In  Corollary  4.2  below, 
we  show  that  this  system  has  a  six-dimensional  solution  space. 

The  combination  l/^^A^^VA0),  which  occurs  in  various  places  in  the  above  for¬ 
mulas,  admits  a  convenient  representation:  i.e.,  combining  (2.2)  and  (2.8)  we  deduce 
that 


(3.19) 


p(°)\(°) 


-V®)  = 


^cos2^ 

—  sin  2if> 

7Icos2^ 

0 

0 

L  71 cos  2^  J 


The  following  supplementary  information  will  be  needed  in  Section  4. 

Lemma  3.1.  Under  the  reflection  9  i— >  — 9 ,  the  functions  in  separation  of  variables 
have  the  parities 


(a)  p(-9)  =  (-l)mm,  (b)  41\-9)  =  (-l)m41\9) 

(c)  v£\-9)  =  (-ir+141\9),  (d)  v$\-9)  =  (- 


(3.20) 
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A2(6) 

MO) 

MO) 

bi(0) 

bo(0) 

di(0) 

do(0) 


i  0  0 
0  |  0 
0  0  i 


+  0(0) 


0  -  — 
U  3 


0 

2  u 

n  _ m2  5 

U  2  6 

r,  _ 4m 

U  TT 


+  0(1) 


0 

_ 4m 

5m^3 _  1 

6  2 


+  0(6- 


(1  + s/a/3)  [010  ]T  +  0(0) 
0_1(l  +  s/-\/3)  [  0  0  -m]T  +  0(  1) 
[  0  1  0  ]T  (exactly) 

6*-1  [  0  1  m  ]T  +  e>(l) 


Table  3.4 

Leading-orders  in  the  expansions  at  9  =  0  of  the  coefficient  matrices  of  (3.17-3.18) 


Proof.  The  reflection  0  i— >  — 6  and  the  rotation  <f>  i— >  (f>  +  7r  are  different  represen¬ 
tations  of  the  same  mapping.  Therefore,  since  p  is  a  scalar 

p( — 0)  cos  m<j>  =  p(0)  cos  m(</>  +  n)  =  (— 1  )mp(0)  cos  m<j), 

which  proves  (3.20a).  Equation  (3.20b)  follows  from  the  same  argument  since  vr 
transforms  as  a  scalar  under  changes  in  the  angular  coordinates.  Rather  than  analyze 
the  parities  of  vg  and  v^,  we  prefer  an  indirect  argument.  Since  V  •  tv1)  is  a  scalar, 
V  •  t/1)  has  parity  (— l)m  under  the  reflection  6  i— >  — 0 ,  and  on  inspecting  (3.7),  we 
deduce  (3.20c,d).  □ 

Incidentally,  although  we  shall  not  need  that  information  below,  we  remark  that 
under  this  reflection  Trr,  Tgg ,  Tg and  have  parity  (— l)m  while  Trg  and  have 
parity  (-l)m+1. 

3.2.  Boundary  conditions  at  the  centerline.  Equations  (3.17,3.18)  have  a 
regular  singular  point  at  0  =  0.  The  leading  orders  of  the  coefficient  matrices  in 
these  equations  are  given  in  Table  3.4.  This  information  may  be  determined  without 
knowing  the  Jenike  solution  explicitly  since,  using  the  fact  that  ip(0)  =  0,  we  deduce 
from  (3.19)  that 


iy(o)  s  T 

„  (0)  =  —=  [  —2  0  1  0  0  if, 

p( °)A(o)  W  V3  L  J 

According  to  the  method  of  Frobenius  [3],  equations  (3.17,3.18)  admit  solutions 
of  the  form 


(3.21) 


vW(G)  =  0VF(9),  p{1)(6)  =  r"1/(6») 
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where  F(6)  and  f{9)  are  analytic  near  9  =  0.  Suppose  the  exponent  v  is  real;  if  1  <  1/ 
such  a  solution  is  continuous,  if  v  <  0  it  is  singular,  and  if  0  <  v  <  1  it  is  continuous 
provided  /( 0)  =  0. 

Proposition  3.2.  There  are  exactly  three  linearly  independent  solutions  of 
(3.17,3.18)  of  the  form  (3.21)  that  are  continuous  at  9  =  0. 

Proof.  Substitution  of  (3.21)  into  (3.17,3.18)  gives  an  indicial  equation  with  roots 

(3.22)  v  =  ±(m  +  1),  ±m,  ir(m  —  1). 

First  suppose  m  >  2.  Since  three  of  the  roots  (3.22)  are  negative,  there  are  at 
most  three  continuous  solutions.  However,  because  the  positive  roots  differ  by  inte¬ 
gers,  these  roots  might  produce  fewer  than  three  continuous  solutions:  the  solutions 
corresponding  to  the  exponents  v  =  m  or  v  =  m  —  1  might  contain  log  terms.  Nev¬ 
ertheless,  using  Maple  we  have  verified  that  this  possibility  does  not  arise — all  three 
values  of  v  produce  continuous  solutions. 

If  to  =  1,  there  are  four  nonnegative  roots  (3.22),  with  zero  being  a  double  root. 
Again  using  Maple  we  have  eliminated  the  various  alternative  possibilities  to  show 
that  there  are,  in  fact,  exactly  three  linearly  independent  continuous  solutions.  See 
Section  5.3  for  more  details  about  the  Maple  code.  □ 

Incidentally,  since  the  roots  of  the  indicial  equation  are  integers,  the  continuous 
solutions  of  the  lemma  are  actually  analytic  near  9  =  0. 

As  noted  above,  we  prove  in  Corollary  4.2  that  equations  (3.17,3.18)  have  a  six¬ 
dimensional  solution  space.  (This  fact  also  emerges  in  the  Maple  computation.)  Thus, 
the  condition  that  solutions  be  regular  at  9  =  0  is  equivalent  to  three  boundary 
conditions.  Therefore,  regularity  at  9  =  0  plus  the  three  boundary  conditions  (2. 5, 2. 6) 
will  provide  a  complete  set  of  boundary  conditions. 

3.3.  Boundary  conditions  at  the  hopper  wall.  We  derive  the  perturbed 
version  of  (2.5)  in  some  detail;  similar  issues  arise  for  (2.6),  and  we  treat  the  latter 
equation  more  succinctly.  The  calculations  are  greatly  simplified  by  the  fact  that  we 
may  neglect  any  quantity  that  is  0(e2).  To  exploit  this  simplification  efficiently,  we 
temporarily  use  the  notation  F  ~  G  to  mean  that  F  =  G  +  0(e2). 

Including  a  prefactor  of  r2  to  remove  all  r-dependence  from  the  equation,  we  may 
rewrite  (2.5)  as 

(3.23)  r2v(r,  9W  +  e  cos  m(f>,  </>)  •  N  =  0. 

Because  of  the  perturbation,  (3.23)  differs  from  (2.9)  in  three  respects: 

-  the  velocity  v  contains  an  additional  term,  v  ~  +  ei/1-*; 

-  the  velocity  is  evaluated  at  a  location  shifted  by  e  cos  m(j>,  and 

-  the  direction  of  the  normal  N  is  changed. 

Regarding  the  first  two  points,  we  observe  that 

r2v(r,9w  +  ecosm<j),<j>)  ~  v^(0w)  +  ecos  m(f>  dgv^°\9w)  +  etrig  m,(f>v^1\9w) 

where  trig  m4>  equals  cos  mcf)  or  sin  m(f>,  depending  on  the  component  of  Regard¬ 
ing  the  third  point,  dQ  is  the  zero  set  of  the  function  9  —  9W  —  ecos m<t>.  Taking  the 
gradient  of  this  function,  we  conclude  that  the  (inward)  normal  is 

TV  ~  [  0  -1  1T  . 

sin  Oyj 
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Modulo  an  0(e2)-error,  N  has  unit  length.  Substituting  the  previous  two  equations 
into  (3.23),  we  deduce  that 


—r2v 


(r,  9w  +  e  cos  rruf>,  <j>)  ■  N  ~  Vg°\ow)  +  e  cos  m<j)  (dgv^  (0W)  +  Vg1)  (0W)J 


esin  m4  (o) 

6  sin  0W  V *  [Uwh 


However,  since  and  vanish  identically  for  Jenike’s  solution,  the  velocity  bound¬ 
ary  condition  for  the  perturbed  problem  reduces  to 


(3.24) 


f\0v,)  =  o. 


We  turn  to  the  stress  boundary  condition  (2.6).  As  regards  the  scalar  tn  in  (2.6), 
we  observe  that,  since  T^1  and  vanish  for  Jenike’s  solution, 


(3.25) 


tn  = 


3 

E 

*0  =  1 


TijNiNj 


p(0) 

lee 


eT 


(i) 


99 


The  vectors  tt  and  vt  in  (2.6)  lie  in  a  two-dimensional  subspace  tangent  to  dtt.  Note 
that  the  unperturbed  tangent  space  is  spanned  by  the  r  and  <f>  coordinate  directions. 
Even  allowing  for  the  perturbation,  the  two  sides  of  (2.6)  will  be  equal  iff  their  r-  and 
(^components  are  equal;  in  symbols,  iff 


TTr 

HwTN 

vr 

TT<j> 

\v\ 

.  V(p 

This  equality  will  hold  iff  (i)  the  two  sides  of  the  equation  are  parallel  vectors  and 
(ii)  the  first  components  of  the  two  sides  are  equal;  again,  in  symbols,  iff 


(3.26) 

(3.27) 


TTr^tj,  —  TT(f,Vr  =  0 
TTr  +  HwTN{Vr/\v\)  =  0. 


and 


Regarding  v,  it  is  clear  that 
(3.28) 


vr 

v<f, 


(0) 

Vr 


J1) 


Regarding  tt  =  r  —  tnN,  we  claim  that 


(3.29) 


TTr 

TT<t> 


T' 


(o) 


T 


(i) 


T, 


r0 

(1) 


Verifying  this  claim  is  straightforward  except  that,  in  analyzing  the  second  component, 
one  must  invoke  the  fact  that  Jenike’s  solution  satisfies  T##  =  .  On  substituting 

(3.28)  and  (3.29)  into  (3.26),  we  obtain  the  equation 

e  ( Tl°ev ^  ~  vl°'>T9^S)  =  0  at  0  =  0W  +  e  cos  m<j). 

The  difference  between  evaluating  this  expression  at  6  =  0W  and  at  the  perturbed  lo¬ 
cation  is  0(e2).  Removing  the  r-dependence  (proportional  to  r)  and  the  ^-dependence 
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(proportional  to  sin  nup)  from  this  equation,  we  obtain  the  first  stress  boundary  con¬ 
dition  for  the  perturbed  problem: 

(3.30)  (t$vW  -  vMf, =  o  ate  =  ew. 

Regarding  (3.27),  we  claim  that 


M- 


Indeed,  it  is  clear  from  (3.28)  that  the  contribution  of  v<f,  to  |u|  is  0(e2),  and  by  (3.24) 
the  contribution  of  vg  to  |u|  is  0(e4).  Thus,  vr/\v\  ~  —1.  Substituting  (3.25)  and 
(3.29)  into  (3.27),  we  obtain  the  condition 

C T- +  eTre'> )  +  Vw (Tgg'1  +  eTgl ) )  ~  0  at  9  =  0W  +  e cos m<t>. 


By  (2.11),  +  HwTgg^  vanishes  at  0  =  6W,  but  at  the  perturbed  location  these 

terms  make  an  0(e) -contribution.  Allowing  for  this  contribution  and  eliminating  the 
r-  and  (^-dependence,  we  derive  the  second  stress  boundary  condition  for  the  perturbed 
problem: 

(3.31)  +  ^wTgg  =  -dg  (f$  +  fiwT^  at  d  =  dw. 


We  have  put  the  inhomogeneous  term,  which  does  not  involve  the  perturbation  T^\ 
on  the  right  side  of  the  equation.  (By  contrast,  (3.30)  and  (3.24)  are  homogeneous.) 

It  is  noteworthy  that  the  perturbed  boundary  conditions  (3.30,  3.31)  resemble 
(2.10,  2.11)  rather  closely. 

4.  Numerical  approximation  of  the  2-point  BVP.  The  coefficients  in  (3.17, 
3.18)  depend  on  the  zeroth-order  solution  discussed  in  Section  2.2.  This  solution  can 
be  found  numerically  without  difficulty,  see  e.g.  [7]  where  a  shooting  method  is  used 
or  [11].  We  will  consider  the  zeroth-order  solution  as  given,  and  we  will  focus  on  the 
corrections  tv1)  and  p^. 

To  simplify  the  notation  before  discretization,  we  set 


=  D«,  *=  40(1)  and  q  =  pM 
at) 


and  rewrite  equations  (3.17,  3.18)  as  a  first-order  system 


'I  0  o' 

'  w'  ' 

\  0  -I  0  1 

w 

'  0 ' 

0  —  A2  b\ 

+ 

0 

-Q 

0 

z 

= 

0 

0  0  0 

<! 

do  df  0 

Q 

0 

where  the  coefficient  matrices  are  the  same  as  above.  The  system  (4.1)  is  completed 
by  the  three  boundary  conditions  (3.24,  3.30,  3.31). 

The  above  system  (4.1)  is  differential-algebraic;  in  the  next  lemma  we  show  it 
has  index  one.  (The  meaning  of  this  term  is  defined  in  the  proof,  or  see  [4].)  The 
approximation  of  solutions  of  the  initial-value  problem  for  such  low-index  DAEs  is 
relatively  well-understood;  see  for  instance  [4]  for  convergence  results.  Moreover,  some 
results  for  the  initial-value  problem  may  be  extended  to  boundary-value  problems,  see 
[5] .  These  considerations  provide  a  theoretical  justification  for  our  using  the  midpoint 
rule  to  solve  (4.1)  numerically. 
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Lemma  4.1.  Assuming  downward  flow,  i.e.,  vr(9)  <  0  for  any  9,  the  first-order 
system  is  differential-algebraic  of  index  one. 

Proof.  We  need  to  show  that  by  differentiating  some  of  the  components  of  (4.1)  at 
most  once,  the  algebraic  character  of  the  system  can  be  eliminated,  leaving  a  purely 
differential  equation.  Let  us  differentiate  only  the  last  component  of  (4.1), 

(4.2)  d,Q  w  +  df  z  =  0. 


The  resulting  system  may  be  written 


o 

o 

'  w’  ' 

linear 

0  —A2  b\ 

z' 

+ 

zeroth-order 

0  df  0 

q' 

terms 

We  claim  the  coefficient  matrix  in  (4.3)  is  nonsingular.  Then,  multiplying  (4.3)  by 
the  inverse  of  this  matrix,  we  obtain  a  purely  differential  equation. 

To  prove  the  claim,  it  suffices  to  show  that 


(4.4) 


+  A(°U2  &! 

dj'  0 


is  nonsingular,  where,  without  changing  invertibility  we  have  inserted  a  factor  of 
in  the  upper  left,  which  simplifies  the  calculation.  Let  us  introduce  the  notation  W 
for  the  column  vector  on  the  RHS  of  (3.19),  so  that  l/^p^A^))^0^  =  sW.  Then 
from  the  definitions  following  (3.17,3.18),  we  have  b\  =  gi  +  DiW.  Similarly,  regarding 
A2,  since  MG±  =  Df ,  we  have  A^°^2  =  DiGi  —  \{DiW){DiW)t .  But 


D{W 


—  sin  2-0 


7Icos2^ 


Hence  the  matrix  (4.4)  equals 


B  = 


\  —  \  sin2  2tp 
cos  2 ip  sin 
0 
0 


*  *  —  ssin2?/> 

*  *  1+^|  cos  2f> 

0  i  0 

1  0  0 


where  *  indicates  elements  that  do  not  affect  the  invertibility  of  B.  It  is  readily 
calculated  that 


det  B  = 


1 

4 


^cos2  2if  + 


75 


As  shown  on  p.43  of  [12],  the  assumption  that  v,0>  <  0  implies  that  |^>(0)|  <  7t/4,  and 
the  claim  follows.  □ 

Corollary  4.2.  The  solution  space  of  (4-1)  has  dimension  six. 

Proof.  The  solution  space  of  (4.3),  which  is  seven  dimensional,  may  be  param¬ 
eterized  by  initial  values  [  w(90)  z(90)  q(90)  ]T.  Since  (4.3)  was  obtained  from 

(4.1)  by  differentiating  (4.2),  we  conclude  that  for  a  solution  of  (4.3), 

d^w{9)  +  dfz(9)  =  0  if  and  only  if  d^  w(9q)  +  dfz(9o)  =  0. 
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Thus  the  solution  space  of  (4.1)  may  be  identified  with  the  set  of  solutions  of  (4.3) 
whose  initial  conditions  satisfy  the  scalar  equation  (4.2).  □ 

The  boundary  value  problem  (4.1),  (3.24,  3.30,  3.31)  is  discretized  using  a  sym¬ 
metric  implicit  Runge-Kutta  method  [2],  [4].  Since  the  solutions  are  expected  to 
behave  smoothly  with  respect  to  9,  the  simplest  of  those  methods,  namely  the  mid¬ 
point  rule,  is  chosen.  In  spite  of  being  only  second  order  accurate,  this  choice  is 
shown  to  be  adequate  below.  The  interval  (0,0^,)  is  divided  in  N  subintervals  of  size 
A 9  =  9W/N ,  defining  a  uniform  mesh  with  nodes  9i  =  i  A 9,  i  =  0, 1,.,. . ,  N.  At  each 
grid  point  9i  there  are  seven  unknowns, 

Ul  =  [w\  w\  w\  z\  z’2  zl3  ql]T. 

Since  there  are  IV +1  grid  points,  there  are  7(iV  +  l)  unknowns  in  total.  The  midpoint 
rule  for  the  ODE  (4.1)  is  applied  on  each  interval  [6i~i,  9i\ ,  i  =  1, . . . ,  N,  leading  to 
7 N  equations  for  the  7 (N  +  1)  unknowns. 

Seven  additional  equations  are  needed  to  close  the  system,  and  these  are  provided 
by  the  boundary  conditions.  At  9  =  9W,  the  three  conditions  (3.24,  3.30,  3.31)  are 
imposed;  and  at  9  =  0,  the  four  numerical  boundary  conditions  listed  in  Table  4.1  are 
imposed.  The  latter  boundary  conditions  may  be  justified  as  follows.  According  to 
(3.21,3.22),  as  9  ->  0, 


W~9V,  q  ~  9l'~1, 

where  v  >  m  —  1.  Thus  u>(0)  =  0  if  m  >  2,  and  q{ 0)  =  0  if  m  >  3.  In  fact,  if  m  =  2, 
direct  calculation  of  the  Frobenius  solution  (3.21)  shows  that  g(0)  =  0  remains  true 
in  this  case  too.  If  to  =  1,  we  refer  to  Lemma  3.1:  by  parity,  w i,  g,  and  Z3  =  dw^/d9 
all  vanish  at  9  =  0.  The  fourth  boundary  condition  in  Table  4.1  follows  from  the  last 
equation  in  (4.1)  in  the  limit  9  — >  0. 


TO.  =  1 

=  0  W2  +  W3  =  0  £3=0  gu  =  0 

TO  >  2 

il’i  =0  w2  =  0  iCg  =  0  q°  =  0 

Table  4.1 

Numerical  boundary  conditions  at  6  =  0. 


The  resulting  7(1V  +  1) 


x  7 (N  +  1)  system  has  the  following  structure 


'  Sj  R\ 

'  U°  ' 

0 

S2  r2 

u 1 

0 

Sn  Rn 

UN-1 

0 

Bo  Bw 

_  UN 

.  Q  . 

The  last  row  of  the  above  system  corresponds  to  the  implementation  of  the  boundary 
conditions;  the  7x7  matrices  B0  and  Bw  contain  the  coefficients  entering  in  the 
formula  from  Table  4.1  and  (3.24,  3.30,  3.31)  respectively,  while  Q  corresponds  to  the 
nonhomogeneous  part  of  the  boundary  condition  (3.31). 
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Fig.  5.1.  Stream  function  showing  secondary  flow  in  a  tilted  hopper  (m  =  1,  9W  =  30°, 
<5  =  30°);  top:  angle  of  wall  friction  =  15°  =  tan  15°),  bottom:  angle  of  wall  friction  =  23.3° 

(fxw  =  tan 23.3°).  By  symmetry,  only  half  of  the  hopper  is  represented. 


5.  Numerical  results. 

5.1.  Secondary  circulation.  We  claim  that,  for  solutions  of  (4.1),  secondary 
circulation — i.e. ,  flow  tangential  to  the  spherical  cap  {r  =  const} — may  be  described 
in  terms  of  the  stream  function 


=  —  sind  sinm</>u>2(6b. 
mr 


In  other  words,  we  must  show  that 

(5.1)  (a)  4°  =  — (b) 

r  suit/  v  r 

Since  1  =  r~2W2(9)  cos  m<j>,  equation  (5.1a)  follows  by  direct  differentiation.  On 
the  other  hand,  since  =  r~2w\(6)  cos  rruj>,  we  have  (dr  +  2r~1)vr  =  0,  so  by  (3.7) 


(dg  +  cot  0)w2(0)  cosmcj)  +  esc  9w^{9)  9^(sin m<t>)  =  0, 


from  which  (5.1b)  follows. 
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Fig.  5.2.  Stream  function  showing  secondary  flow  in  an  “elliptical”  ( m  =  2,  0W  =  30°,  6  = 
30°);  top:  angle  of  wall  friction  =  15°  (//„•  =  tan  15°),  bottom:  angle  of  wall  friction  =  23.3° 
=  tan  23.3°). 


Figures  5.1  and  5.2  show  plots  of  the  level  lines  of  'F,  which  equal  the  projection 
of  the  streamlines  onto  a  spherical  cap  {r  =  const}.  Figure  5.1  corresponds  to  a  tilted 
hopper  (to  =  1),  while  Figure  5.2  corresponds  to  an  “elliptical”  hopper  (to  =  2). 
In  both  cases,  two  different  values  of  the  wall  friction  are  presented,  while  the  other 
parameters  are  held  constant.  The  grains  do  not  move  along  radial  lines  but  follow 
more  complicated  and  fully  three-dimensional  trajectories.  One  can  observe  from 
Figures  5.1  and  5.2  that  the  sign  of  the  main  circulation  changes  when  /iw  increases. 
In  the  case  to  =  1,  this  is  even  accompanied  by  a  change  in  the  topology  of  the 
flow — i.e. ,  the  addition  of  a  vortex. 

As  varies  between  the  values  shown  in  the  Figures,  the  flows  do  not  undergo 
a  smooth  transition  from  one  case  to  the  other.  This  is  demonstrated  in  Figure  5.3 
which  shows,  for  values  of  to  from  1  to  4,  the  azimuthal  velocity  v^\0w)  on  the  wall 

as  a  function  of  the  wall  friction  fj,w.  For  each  value  of  to,  1  suffers  a  “1/x-blowup” 
as  /iw  passes  through  a  critical  value. 

As  fiw  crosses  the  critical  values  in  Figure  5.3,  the  direction  of  the  circulation 
changes  in  a  singular  way.  It  turns  out  that  there  are  also  additional  critical  parameter 
values  for  which  the  sign  of  the  circulation  changes  smoothly,  passing  continuously 
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wall  friction  u 

Fig.  5.3.  Dependence  of  the  resonance  on  the  geometry  of  the  domain  through  the  coefficient 
m:  blowup  of  (6W)  as  a  function  of  the  wall  friction  fxw  (internal  friction  5  =  30°,  half  opening 
angle  6W  =  30°). 

through  zero.  For  the  case  m  =  1  and  S  =  30°,  curves  of  9W,  pw  along  which  the 
circulation  changes  sign  by  either  mechanism  are  shown  in  Figure  5.4.  Note  that 
the  smooth-transition  curve  does  not  depend  on  the  value  of  to,  but  is  a  property 
of  the  radial  solution  itself.  Specifically,  the  circulation  vanishes  when  the  boundary 
condition  for  the  correction  terms  (3.31)  is  homogeneous,  i.e. ,  d$T$  +  pwdgT^  =  0 
at  9  =  9W.  The  range  of  9W  in  Figure  5.4  is  limited  by  the  mass- flow  limit — exceeding 
this  limit  leads  to  flows  with  rigid  regions,  to  which  the  present  model  does  not  apply. 
The  range  of  fiw  is  limited  by  the  condition  that  /iw  <  sin  S  =  1/2;  here  the  upper 
bound  corresponds  to  a  fully  rough  wall  [8]. 

Figure  5.5  offers  a  three-dimensional  view  of  which  combinations  of  the  parameters 
S,  pw  and  9W  lead  to  blowup.  Another  surface  of  critical  values  corresponds  to  the 
above  mentioned  smooth  transitions.  In  Figure  5.5,  the  curve  of  intersection  between 
those  two  critical  surfaces  is  also  represented. 

5.2.  Relation  to  bifurcation  theory.  While  the  blowup  in  the  correction  so¬ 
lution  could  not  have  been  anticipated,  its  cause  is  easily  understood  a  posteriori.  The 
two-point  boundary  problem  for  T W  has  inhomogeneous  boundary  conditions, 
but  at  the  singular  point,  the  problem  with  the  corresponding  homogeneous  boundary 
condition 

(5.2)  +  pwT$  =0  at  9  =  9W 

has  a  nonzero  solution.  This  behavior  is  analogous  to  that  of  a  forced  harmonic 
oscillator 


x  +  u>qX  =  A  sin  cot. 

The  steady  state  oscillations  have  amplitude  A /(to g  —  w2)  which  diverges  as  ui  passes 
through  the  natural  frequency  uq  of  the  oscillator.  Unlike  the  above  problem,  in  our 
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Fig.  5.4.  Critical  values  leading  to  sign  changes  of  the  circulation  (internal  friction  S  =  30°, 
m  =  1). 
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0.28 
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Fig.  5.5.  Locus  of  parameters  for  which  blowup  of  v occurs  (m  =  1.)  The  dark  curve 
represents  the  intersection  of  this  surface  with  the  locus  of  parameters  for  which  the  flow  changes 
direction  smoothly.  Compare  with  5.4. 


model  the  inhomogeneity  is  in  a  boundary  condition  rather  than  the  equation. 

These  remarks  suggest  a  connection  with  bifurcation  theory.  Consider  equations 
(2. 1-2. 3)  on  the  conical  domain  {0  <  9  <  9W}  subject  to  boundary  conditions  (2.9- 
2.11).  Let  us  write  W  =  [v.  T,  A]T  as  a  multicomponent  unknown,  and  from  the  three 
parameters  <5,  /iVJ ,  0W  let  us  consider  fxw  as  a  distinguished  (bifurcation)  parameter. 
We  rewrite  (2. 1-2. 3),  (2.9-2.11)  symbolically  as 


(5.3) 


=  0. 
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Now  for  any  fiw,  the  Jenike  solution  provides  one  solution  of  (5.3).  We  investigate 
other  possible  solutions  of  this  equation  with  the  implicit  function  theorem.  Let  L  be 
the  linearization  of  $  with  respect  to  its  first  argument:  i.e., 

LW=  lim  +  hW,fiw)  -$(W(0),^)^- 

If  L  is  invertible,  then  the  Jenike  solution  is  the  unique  solution  of  (5.3)  in  some 
neighborhood  of  W^°\  Thus,  to  detect  possible  bifurcation,  we  look  for  values  of  fxw 
such  that  there  is  a  nonzero  solution  of  the  equation 

(5.4)  LW  =  0. 

Essentially,  we  have  already  calculated  L.  In  linearizing  the  PDEs  (2. 1-2. 3)  we 
obtain  (3. 3-3. 5).  In  linearizing  boundary  conditions,  (2.9)  becomes  (3.24);  (2.10)  be¬ 
comes  (3.30);  and,  because  the  location  of  the  boundary  is  not  moved,  (2.11)  becomes 
not  (3.31)  but  (5.2).  By  the  above  analysis,  (5.4)  has  a  nonzero  solution  precisely 
when  fiw  equals  a  critical  value  at  which  the  correction  solution  blows  up. 

In  light  of  bifurcation  theory  [6],  these  observations  lead  us  to  make  the  follow¬ 
ing  conjecture:  Even  in  a  conical  domain,  there  are  solutions  of  (5.3)  with  nonzero 
circulation.  These  bifurcate  from  the  Jenike  solution  where  fiw  equals  a  critical  value. 
We  shall  explore  this  conjecture  in  a  future  publication. 

5.3.  Checks  on  the  computation.  For  comparison  with  the  above  numerical 
solution,  the  method  of  Frobenius  was  applied  directly  to  the  system  (3.17,3.18)  using 
Maple.  Given  Jenike ’s  radial  field,  a  linear  system  for  the  coefficients  of  the  series 
solution  is  readily  formed  and  solved,  yielding  a  solution  with  three  free  parameters, 
corresponding  to  the  three  linearly  independent  solutions  in  Proposition  3.2.  Subse¬ 
quently,  the  three  boundary  conditions  (3.24,3.30,3.31)  provide  the  needed  relations 
to  determine  the  solution  to  the  full  boundary  value  problem. 

Two  methods  of  obtaining  the  radial  field  were  employed.  Under  the  assumption 
that  9 ^  and  nw/9w  are  both  small  and  of  the  same  order,  a  series  representation  of  the 
Jenike  field  was  computed  within  Maple  itself.  Under  the  less  restrictive  assumption 
that  only  9W  be  small  (say  10°),  numerical  solutions  were  computed  in  MATLAB, 
fitted  to  polynomials,  and  then  imported  into  Maple.  In  both  cases,  the  resulting 
polynomials  were  then  used  to  compute  the  first  order  correction.  The  corrections 


x  10‘3 


0  0 


Fig.  5.6.  Comparison  of  from  the  purely  numerical  method  of  Section  4  and  from  the 
Frobenius  method  of  Section  5.3.  (Using  m  =  1,  0W  =  10°,  S  =  30°,  and  [iw  =  0.3.) 
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to  the  stress  and  velocity  obtained  through  this  symbolic  approach  agree  extremely 
well  with  the  results  of  the  purely  numerical  method  of  Sections  4  and  5:  for  the 
representative  values  6W  =  10°,  6  =  30°,  and  fiw  =  0.3  the  corrections  obtained  by 
the  two  different  methods  have  a  relative  difference  of  less  than  1%.  Furthermore, 
with  6W  =  10°  and  5  =  30°  the  corrections  obtained  via  Maple  exhibited  the  familiar 
“l/rc-blowup”  near  /j,w  =  0.4605.  For  comparison,  in  the  purely  numerical  approach, 
at  the  same  values  of  S  and  0W ,  blowup  occurred  for  fiw  =  0.4610. 
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